Phys. Rev. A 74, 052308 (2006) 






% 



quant -ph/0604126 



Genetic Algorithm Optimization of Entanglement 

Jorge C Navarro-Munoz* , H.C. Rosu^ and R. Lopez-Sandoval"'" 
Is I ^ Potosinian Institute of Science and Technology, 

^ ■ Apartado Postal 3-74 Tangamanga, 78231 San Luis Potosi, Mexico 

^ ■ (Received 19 April 2006; published 7 November 2006) 

^ ; (Dated: Nov. 13, 2006) 

m . 

Abstract 

\^J , We present an application of a genetic algorithmic computational method to the optimization 

jH \ of the concurrence measure of entanglement for the cases of one dimensional chains, as well as 
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I. INTRODUCTION 



Entanglement has been described as an important characteristic for quantum informa- 



tion and quantum computation [1|. In the following, the quantity of interest will be the 
concurrence measure of entanglement introduced by Wootters [2!]. It will be studied for 
various kinds of condensed-matter lattice systems and the main focus will be on calculating 
optimized entangled states using genetic algorithms [3] . To the best of our knowledge, at the 
present time there exists only a paper by Prashant on the application of genetic algorithms 



a 



to evolving quantum circuits |4|. 

The motivation for our work resides in recent studies about maximum nearest-neighbor 
entanglement QQ. In these cases, a N-,.Ut Hng in a t..ans.ationa.,v n^anant quantun. 

state has been considered. Under certain conditions, O'Connor and Wooters [5j have found 
formulas to obtain the maximum possible nearest-neighbor entanglement. Moreover, they 
have compared this quantity with the entanglement produced off an antiferromagnetic state 
of a ring with an even number of spin-| particles. Also, there have been studies of con- 
currence for nearest-neighbors in finite clusters with the purpose to see its behavior in two 







dimensions. In particular, this was carried on for square, triangular and Kagome lattices [7|. 
Further studies focus on systems with higher order entanglement, that is, when subsystem 
A is bigger than two qubits 8|. 

In the lattice systems of this research, qubits are represented by sites in a lattice or a 
chain. The two computational basis |0) and |1) are represented by occupied and empty 
sites, respectively. Using this representation, the concurrence can be calculated for different 
fillings. This approach might be useful for physical experiments involving electron control 
(e.g., quantum dots [9|). 

The electronic system will be described by a tight-binding Hamiltonian of the form 

H = Y^ Eihi + ^ tijc\cj , (1) 

where, for simplicity, we will consider spinless electrons. In ([T]) q [cj] is the usual creation 
(annihilation) operator of a spinless electron at site i, whereas hi = cjci is the number 
operator, and tij is the hopping integral between nearest-neighbors (NN) and next-nearest- 
neighbors (NNN) sites i and j. Si is the on-site energy for atom i. For simplicity, we work 
with the same kind of atoms and we take e^ = 0. There exist physical systems which can 



be modeled using a similar approach as here, most notably the polyacetylene systems [IC 
in the static approximation (without phonons). 

Notice that we are using a one-body Hamiltonian. For this case the ground state wave 
function can be obtained simply by filling the lowest one-body eigenstates. The correspond- 
ing eigenvalues are obtained by diagonalizing the matrix Hamiltonian as expressed in the 
one-body basis. In this case, one can always find a linear transformation of c] (q) leading 
to a diagonal form of the Hamiltonian. 

Specifically, our tight-binding Hamiltonian H = — ^ . . tijC^Cj can be written as H = 
J2k^kCkCk, where q = J^h'^ik^^l ^^^ ^i — J2k'^ik^k- For periodic systems with only NN 
hopping integrals (tjj = t) [11], the Hamiltonian H can be diagonalized through the Fourier 
transformations q = -^ X]fc=i^^P(^2^^^/^)^fc ^^^ ^^^ complex-conjugated counterpart, 
therefore H takes the form 
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H = -2t'^cos{27Tk/N)cl 
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On the other hand, for a nonperiodic system with arbitrary hopping integrals it is not 
possible to get an analytic diagonalization procedure. For the simplest nonanalytic case of 
a system of A^ atomic sites and nearest-neighbor hopping integrals, one should diagonalize 
a Hamiltonian matrix of tridiagonal form: 
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We have used the LAPack subroutines 12| to diagonalize this type of Hamiltonian matrices 
through the QR algorithm as well as more complicated forms resulting from next-nearest 
interactions, where the Householder reduction to tridiagonal forms is applied first. The 
eigenvectors for K fermions can be obtained in a direct way using the preceding equations 
and are given by the tensorial product of the one-body eigenvectors. 
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Thus, it is clear that for the ground state we require the lowest levels to be occupied. 

Moreover, through a Jordan- Wigner transformation the spinless fermion Hamiltonian can 
be rewritten as a XX spin-i chain. In the spinless fermion case, each lattice site is either 
occupied or free, whereas in the spin polarization case each lattice site can have the spin 
up or down. It is well known that oxides and fluorides of transition metals, e.g., MnO, NiO 
and MnF2, FeF2, C0F2, respectively, are described by such simple spin Hamiltonians J13l |. 

Wootters' formula that we used to calculate the concurrence is obtained in the Appendix 
and is given by 

C = max{0, 2|p23| - S^/pHpil} = 2 max{0, |p23| - VpuPm} ■ (5) 

The concurrence calculations in Sec. II are always a sum over all the pairs of sites divided 
to the total number of sites. Section III contains some conclusions and the Appendix is 
devoted to the mathematics of the concurrence formula. 

II. OPTIMIZING ENTANGLEMENT USING GENETIC ALGORITHMS 

We pass now to the main goal of the paper which is the maximization of entanglement 
using genetic algorithms. Specifically, we will consider the ground state of a spinless system 
modeled by the tight binding Hamiltonian given in Eq. ([T]). 

We recall that genetic algorithms (GAs) were invented by John Henry Holland in the 
1960s and were developed by him and his students and colleagues at the University of 
Michigan in the 1960s and the 1970s. Holland's goal was not to design algorithms to solve 
specific problems, but rather to formally study the phenomenon of adaptation as it occurs 
in nature and to develop ways in which the mechanisms of natural adaptation might be 
imported into computer systems. 

Much alike nature, Holland's GA is a method for moving from one population of "chro- 
mosomes" (e.g., strings of characters or numbers) to a new population by using a kind of 
"natural selection" together with the genetics-inspired operators of crossover, mutation, and 
inversion (this last operator is rarely used nowadays). 

The genetic pseudoalgorithm employed by us here goes as follows: 

1. Read input parameters including type of lattice, sites in the system, number of gener- 
ations, crossover probability, mutation probabihty, etc. 



2. Build a table with indices of the nearest neighbors of each site. A table including also 
next-nearest neighbors can be built as well. 

3. Using the neighbor table, identify the specific places in the Hamiltonian matrix where 
"bonds" occur. Each place represents a valid tij entries and will be stored in a special 
array. This array will be considered hereafter as a chromosome. 

4. Allocate two arrays, "generationO" and "generationl" composed of chromosomes. 

5. Construct an additional chromosome called "best" with initial random numbers be- 
tween (0, 5). 

6. For a given range of filling repeat: 

• Initialize "generationO" with random values in the range (0,5). 

• Make the first chromosome of "generationO" equal to "best" . 

• For a given number of generations repeat: 

— Decode each chromosome in "generationO" into a Hamiltonian matrix, diag- 
onalize it and calculate the total concurrence between all nearest neighbors 
of the system. In other words, calculate fitness for each individual in "gen- 
erationO" . 

— Make "best" equal to the chromosome with highest value of fitness in "gen- 
erationO" 

— Print the value of the average fitness of the population of "generationO" and 
fitness of "best" in output files. 

— Apply the selection operator: Use crossover and mutation operators on chro- 
mosomes in "generationO" to create new chromosomes into "generationl" . 

— Make "generationO" equal to "generationl". 

— Make the first chromosome in "generationO" equal to "best". 

• Find the chromosome with the maximum fitness. Print its fitness value in an 
output file. 

To make the calculations tractable the biggest system that we considered was of 49 sites 
with 800 generations for which the optimization procedure has taken about three days. In all 



calculations we have worked with the probability of crossover pc = 0.70 and the probability 
of mutation pm = 0.002. 

A. One-dimensional chains 

We begin the analysis of entanglement maximization using genetic algorithms with the 
simplest case of small lineal chains with and without periodic boundary conditions. In 
Figs. [1] and [2] we present results of concurrence as a function of percentage filling for two 
chains with 24 and 44 sites, respectively. Besides nearest-neighbor interactions, we have 
also considered interactions with both nearest neighbors and next-nearest neighbors in the 
Hamiltonian. The population size remained at 400 individuals and the generations were 
kept at 250. Later on, the role of the number of generations will become apparent. 

From the figures, it can be noticed that in the case of nearest-neighbor interactions with 
and without periodic boundary conditions, concurrence as function of filling is smoother 
than in the cases where next nearest neighbors are also considered. This can be due to 
a larger size in the chromosomes in the latter case and a greater number of generations 
are necessary to obtain a similar behavior than its only- nearest-neighbors counterpart. We 
can only conclude that a greater number of generations and possibly a greater size in the 
population is necessary to overcome these oscillations. 

Also, notice that cases including next nearest neighbors cannot yield lower results than 
the only nearest neighbors case. This is because the chromosomes from the former case 
contain the chromosomes of the latter (i.e., the NN case is a subset of the NNN case), which 
offers the possibility to explore a wider spectrum of solutions. In the case where this extra 
space yielded only lower results, the best chromosomes would be those of the NN space. 
This phenomenon can be most clearly noticed near half filling. Once again, this behavior is 
a consequence of a greater number of generations. 

At this point, it is important to remember that there are various parameters responsible 
for a larger chromosome in this kind of system. These parameters are the size of the system, 
the periodic boundary conditions and bringing next-nearest-neighbors interactions into play. 
A larger chromosome would allow an exploration of a wider solution space but on the other 
hand it is expected to decrease the convergence time. 

We have already mentioned the possibility of a greater number of generations affecting 



directly the smoothness of the concurrence. We addressed this question by running two cases 
depicted in Figs. [3] and HJ where the former does not consider periodic boundary conditions 
while the latter does. 

In both cases we have set a 44-site chain with a population size of 400. Only the inter- 
actions between nearest neighbors were taken into account. 

Both figures confirm our early supposition about increasing the number of generations 
since the correlation functions look increasingly smoother. Notice, however, that certain 
roughness still remains. Some possible solutions consist of increasing the size of the pop- 
ulation, dynamically change the mutation probability (when variation between individuals 
begins to narrow) and raising the number of generations further more. Different selection 
methods could also be considered because an inefficient parent selection could lead to slow 
evolution of the system. Even though it is clear that individuals with better fitness are 
obtained, notice how in Fig. [3] the best chromosome near 0.05 filling was obtained with 600 
generations despite having cases with up to three times more generations. 

In Fig. [5] (for its more detailed structure up to 2000 generations see Fig. [6]) we follow 
the evolution (optimization) of concurrence for each filling in a 44-site chain with periodic 
boundary conditions and interactions only between nearest neighbors. The population size 
for this calculation was 350 and the number of generations was 500 per filling. The average 
fitness per population is compared with the fitness from the best chromosome in the popula- 
tion. Notice how the population always follows closely the evolution of the best chromosome. 
Transitions between different fillings are readily noticed through a drop in average fitness. 
A very remarkable feature is that the best chromosome for a certain filling ranks high for 
the next filling but is not the highest. In other words, there are different best chromosomes 
for different band fillings. We also notice that the bottom dots correspond to the average 
concurrence for randomly disordered populations and that the average concurrences for the 
subsequent optimized GA populations are always better than the disordered cases. 

Another remarkable characteristic about Fig. [5] is its symmetry around half-band filling. 
This property is due to the fact that this is a bipartite lattice and consequently its physical 
properties are symmetric because of an electron-hole transformation. The fact that the 
results presented show this property reassures the validity of our calculations. 



B. Two-dimensional systems 

We will study now the optimization of concurrence in two-dimensional systems modeled 
by means of the same tight-binding Hamiltonian. 

1. Square lattices 

In Figs. [7] and [8] we display the concurrence as a function of band filling for a 7x7 square 
lattice. 

Figure [7]presents a comparison between systems using nearest-neighbor and next-nearest- 
neighbor interactions as well as periodic and open boundary conditions. The number of 
generations for these cases has been chosen 600 and the population size 350. It is worth 
mentioning that in general the cases with interactions only between nearest-neighbors rank 
slightly higher in its concurrence value. This is a somewhat unexpected result that may 
be attributed to various factors including selection methods and number of generations. 
Possible reasons for this behavior were addressed in the preceding section. 

To study the effect of the number of generations on the optimized value of concurrence 
and the smoothness of the curve, we present calculations for four different cases in Fig. [8l In 
these cases, population size was kept at 400. It is clear that by raising this number we are 
able to obtain better optimized solutions and the concurrence curve tends to be smoother. 

2. Triangular lattices 

Finally, we have made calculations for not bipartite lattices in order to study the effect 
of frustration on concurrence. It has already been mentioned that these kinds of lattices are 
not symmetric under an electron-hole transformation. This is the reason why their physical 
properties differ completely between lower and upper sections of band filling. 

As a particular case of a nonbipartite lattice, we have considered a triangular lattice with 
49 sites. Results of our calculations are shown in Figs. [3 and dUl 

In Fig. Oa population size of 350 has been used and the system has been allowed to go up 
to 600 generations. As in the preceding sections, this case includes the interaction between 
nearest and next-nearest neighbors, as well as open and periodic boundary conditions. Once 
again, we find better optimizations for nearest-neighbor interactions. It is important to 
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remember that we are dealing with a more complex chromosome, as sites in this kind of 
lattice have a greater number of neighbors than one dimensional systems. This is also a 
cause for a lower time in convergence as the solution space increases considerably. 

The effect of the number of generations can be examined in Fig. [TUl These results 
demonstrate the slow convergence when calculating this kind of system. Notice that, al- 
though oscillations decrease and better individuals are found, efficiency narrows between 
the cases with 600 and 800 generations. 

III. CONCLUSIONS 

We have implemented computational techniques -more specifically genetic algorithms- to 
optimize entanglement in systems modeled according to a tight-binding Hamiltonian. The 
qubits in all these studies have been described as sites in the system and the computational 
basis as occupied or empty sites. 

Our application of genetic algorithms has proven to be valuable, since we obtained 
configurations which yield better results for concurrence in the randomly disordered one- 
dimensional case as discussed in Sec. II A (in fact, we have partial results confirming this 
statement for two-dimensional cases as well). Moreover, the GA optimization provided bet- 
ter results even with respect to the ordered cases as can be noticed in Figs. [TT], [T21 and 

m 

We finally mention that the optimal Hamiltonian for the one-dimensional case with peri- 
odic boundary conditions corresponds to dimerized chain for which the coupling coefficients 
take alternate high-low values in the half filling region. This structural phase transition 
has been found in the polyacetylene chains and suggests that optimal entanglement can be 
obtained in the dimerized phase of the polyacetylene. We have also some evidence of this 
behavior in two-dimensional systems. These results will be reported elsewhere. 

Quantum computation and quantum information are still a long way to go. Nevertheless, 
these areas seems to be a logical and necessary step in tomorrow's technological world. In 
this scenario, quantum entanglement will play a critical role, and our work attempts to be 
another step towards better understanding it. 



APPENDIX: THE CONCURRENCE 

To calculate the concurrence we need first the density matrix pA for qubits i, j, which is 
the trace over system B of all the possible states \'iPab){''Pab\- The general state function 
for this system is 

I^As) = J^a„|^A)|^B) (6) 

n 

where, for a system of A^ sites, n goes through all the 2^ possible combinations in the 
computational basis (e.g., |00 . . . 00) -^ 1 11 ... 11)). Subsystem A is comprised of the two 
qubits of interest in the sites i, j (i.e., iV'j) ® |'?/'j), see also Fig. [H]). For a specific system of 
A^ sites, there are Ni occupied and N — Ni not occupied sites. Our two-qubit subsystem A 
has, naturally, four possible states -namely |00), |01), |10), and |11)- therefore Eq. ^ can 
be decomposed in the following manner: 

I^ab) = $^a„|00) ® 1^^) + 5^&o|01) ® 1^1,) + 5^cp|10) ® \rB) + $^rf,|ll) ® 1^^). (7) 

m o p q 

In this equation, the sums run for all the possible combinations in the \'iPb) space such that 
the number Ni of occupied sites is preserved. For example, if |V^a) = |01), system B is left 
with A^^i — 1 occupied sites. 

To obtain the reduced density matrix it is necessary to perform the trace over system B, 

2N-2 

PA=Y1 ((^1 ® {^b\)\^ab){^ab\{\I) ® 1^^)). (8) 

1=1 

It is clear that applying this operation will not eliminate those terms whose elements in the 

B subsystem in \4'ab){'^ab\ have the same number of occupied sites. The terms in system 

A that are left after the trace operation are of the kind |00)(00|, |01)(01|, |01)(10|, |10)(01|, 

|10)(10|, and |11)(11|. 

The |00)(00| element is spared after the trace operator because its \iPb) elements contain 
the same quantity of occupied sites (i.e., Ni sites). This is a similar case with the |11)(11| 
elements where the \ip%) wave functions contain A^^i — 2 occupied sites. 

In the case of |01)(01|, |01)(10|, |10)(01|, and |10)(10| elements, notice how their \iIj°/) 
wave functions have the same number of occupied sites {Ni — 1). 
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Finally, the elements in the reduced density matrix are 

/pii \ 



Pa 



(9) 



P22 P23 

P32 P33 

V P44/ 

For Pa to be a valid density matrix, it must be Hermitic [pA = Pa) ^'^^ its trace be equal 
to 1. This means that P32 = P23 and pn + P22 + P33 + Paa = 1 so it is necessary to calculate 
only four elements of the matrix. 

In order to calculate each of these reduced density matrix elements, the second quanti- 
zation approach will be used. 

The first element of the matrix, pn can be realized as follows 

Pii = (^ab| (1 - fii) (1 - hj) \i)AB) , (10) 

where the operator fij projects on all the elements of the type |xl) ® \iPb) and after applying 
(1 — fij) we end up with all the elements that do not occupy the site j (i.e., |a;0) (X> iV^s))- A 
similar approach follows (1 — Uj) and after applying the bra operation we are left only with 
the coefficients of all the |00) \iPb) states. 

The other elements are obtained likewise with the following operators: P22 = ('^ab|(1 ~ 
hi)hj\2ljAB), P33 = {i^AB\ni{l - nj)\tpAB), P44 = {i^AB\nihj\iJAB) , and P23 = H^abIcjcIUjab)- 

In the matrix element P23, c] maintains only those states of the form |0a;) \^b) trans- 
forming them into |la;) (g) \iPb)- Out of this set of states, Cj deletes all states of the type 
|xO) \iPb) and we end up with states |10) (g) \iPb)- 

It is very easy to show that the pA elements can be calculated as average quantities of 
the complete ground-state wave function. For example, 

Pll = {i^AB\i^AB)-{i^AB\ni\i^AB)-{i^AB\nj\i^AB) + {i^AB\nihj\tpAB) = I- (fli) - (hj) + {fliflj) . 

(11) 

The other elements are obtained similarly, P22 = {fij) — {hihj), P33 = {hi) — {hihj), P44 = 
{hihj), P23 = {cjc}). 

Concurrence is an entanglement monotone in its own right (i.e., positive or zero for any 
density matrix p; for factorizable states and 1 for the Bell states). A simple formula for 
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the concurrence has been worked out by Wootters in 1998 



C{p) = max{0, Ai - A2 - A3 - A4}, 



(12) 



where the A coefficients are the square roots of the eigenvalues of the non-Hermit ian matrix 
PaPa in decreasing order. The formula applies for the density matrix of the subsystem 
with the pair of qubits [pa ='^t^b{p)]- The density matrix p is defined through a spin Rip 
transformation expressed in terms of the imaginary Pauli matrix ay as follows 



Pa 



( 



(Jn 



s)pa( 



cr„ 



ay). 



(13) 



This leads to 



/ -1 \ / p*i 



pa 



10 

10 

\-l / 



\ 



P22 P23 
P32 P33 



/ -l\ 
10 
10 



yo op:4/V-iooo/ 



Now we are able to construct the non-Hermitian matrix PaPa 



( 



PaPa 



\ 



Piipli 

P22P33 + P23P23 P22P32 + P23P22 

P32P33 + P33P23 P32P32 + P33P22 




I pIa ^ \ 

P33 P32 

p;3 P22 

V Q p\j 



\ 



(14) 



(15) 







P11P44 / 

However, pA is indeed Hermitian so the following relationships are taken into account: pn = 
Pii) P22 = P22; P32 = P235 P33 — P33; ^^d P44 — PIa- Therefore, the matrix PaPa has the form 



/ 



PaPa 



\ 



P11P44 










P22P33 + IP23I 

2P33P23 







2p22P23 



\ 



(16) 






P22P33+ |P23|^ 

P11P44 / 

In a block diagonal matrix, the eigenvalues are simply the eigenvalues of individual blocks, 
so two eigenvalues are readily available. The other two are obtained from the following 
determinant 

P22P33 + IP23I " ^ 2p22P23 



2P33P23 



P22P33 + IP23I 



A 



(17) 
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This gives 

A = P22P33 + |p23p =F 2VP22P33IP23I • (18) 

Thus, the four possible values of the A coefficients are as follows: 

Aa = (VP22P33 - |P23|)^ Ab = (VP22P33+ |P23|)^ Ac = P11P44, ^d = PllPH- (19) 

Finally, the square roots of these lambda coefficients are directly employed in ( fT2l) . 

Noticing that A^ is the largest eigenvalue, the latter results leads immediately to Wootters' 
formula Eq. (2). 
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FIG. 1: (Color online) Optimized concurrence as a function of band filling for a linear 

chain with 24 sites. 250 generations; population = 400. Note that in all figures the 

quantities plotted are dimensionless and the values of the parameters pc and pm are fixed 

at 0.70 and 0.002, respectively. 
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FIG. 2: (Color online) Optimized concurrence as a function of band filling for a linear 
chain with 44 sites. 250 generations; population = 400. 
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FIG. 3: (Color online) Optimized concurrence for different number of generations in a 
one-dimensional 44-site lattice. Open boundary conditions, population size = 400. 
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FIG. 4: (Color online) Optimized concurrence for different number of generations in a 
one-dimensional 44-site lattice using periodic boundary conditions, population size = 400. 
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FIG. 5: (Color online) Best and average fitness for each generation in a one- dimensional 
44-site lattice subjected to periodic boundary conditions, nearest-neighbor interactions, for 

500 generations, population size of 350. 
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FIG. 6: (Color online) Detail of Fig. [5] for the first 2000 generations. 
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FIG. 7: (Color online) Optimized concurrence for a square lattice comparing 
nearest-neighbor interactions, next-nearest-neiglibor interactions and boundary conditions. 

600 generations; population size = 350. 
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FIG. 8: (Color online) Optimized concurrence for a square lattice comparing number of 
generations for nearest-neighbor interactions, periodic boundary conditions, population 

size = 400. 
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FIG. 9: (Color online) Optimized concurrence for a triangular lattice comparing 
nearest-neighbor interactions, next nearest neighbor interactions and boundary conditions. 

600 generations; population size = 350. 
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FIG. 10: (Color online) Optimized concurrence for a triangular lattice comparing number 
of generations for nearest-neighbor interactions, periodic boundary conditions, population 

size = 400. 
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FIG. 11: (Color online) Concurrence for the optimized and ordered cases in a 49 sites 
square lattice using nearest-neighbor interactions and periodic boundary conditions 800 

generations; population size = 400. 
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FIG. 12: (Color online) Concurrence for the optimized and ordered cases in a 49 sites 
triangular lattice under the same conditions as in the previous figure. 
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FIG. 13: (Color online) Comparison between the optimized and the ordered case in a 
one-dimensional system of 49 sites. The calculation parameters were identical to those of 

Figs, [m and [HI 
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FIG. 14: (Color online) Schematic illustration of the partition of the system of interest into 
two subsystems for the calculation of the concurrence between sites i and j. 
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